###################
## Load packages ##
###################
library(tidyverse)
library(ggplot2)
library(ggthemes)
library(ggpubr)
library(readxl)
library(patchwork)
library(GGally)
library(reshape2)
library(countrycode)
library(fect)




##################################################
## Set working directory / create output folder ##
##################################################
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
dir.create("output")
dir.create("output/results")




###############
## Load data ##
###############
main_dat <- readRDS("data.RDS")
coder_dat <- readRDS("data_coder.RDS")

vdem <- readRDS("vdem13.RDS")
acd <- readRDS("acd.RDS")
ciri <- readRDS("ciri.RDS")
cv <- readRDS("cv.RDS")
pts <- readRDS("pts.RDS")




##############
## Figure 1 ##
##############
cp_ed <- main_dat %>%
  select(year, epsm_cent, v2xed_ed_cent) %>%
  na.omit %>%
  group_by(year) %>%
  summarize_all(mean, na.rm=T) %>%
  rename(EPSM=epsm_cent,
         `V-Indoc`=v2xed_ed_cent) %>%
  pivot_longer(-year) %>%
  mutate(title="Global Mean: Education Centralization") %>%
  ggplot(., aes(x=year, y=value, group=name, color=name, linetype=name)) +
  theme(plot.title=element_text(size=12, hjust=0.5),
        plot.margin=unit(c(0.5, 0.75, 0.5, 0.75), "cm"),
        axis.title.y=element_blank(),
        axis.text.y=element_text(size=11, color="black"),
        axis.title.x=element_blank(),
        axis.text.x=element_text(size=11, color="black"),
        legend.title=element_blank(),
        legend.key=element_rect(colour=NA, fill=NA),
        legend.key.width=unit(1, "cm"),
        legend.text=element_text(size=11),
        legend.position="bottom",
        legend.margin=margin(0,0,0,0),
        legend.box.margin=margin(-5,-5,-5,-5),
        panel.border=element_rect(colour="black", fill=NA, linewidth=1),
        panel.background=element_rect(fill="white"),
        panel.grid.major=element_line(color="grey", linewidth=0.5),
        strip.background=element_rect(color="black", fill=NA),
        strip.text=element_text(size=11, hjust=0.5)) +
  geom_line(linewidth=0.75) +
  scale_linetype_manual(values=c("solid", "dashed")) +
  scale_color_manual(values=c("blue4", "springgreen3")) +
  scale_x_continuous(breaks=seq(1950, 2020, 20),
                     limits=c(1945, 2020),
                     expand=c(0,0)) +
  scale_y_continuous(breaks=seq(0,1,0.25),
                     limits=c(0,1),
                     expand=c(0,0)) +
  facet_grid(.~title)


cp_dem <- vdem %>%
  select(year, v2x_regime, e_boix_regime) %>%
  mutate(v2x_regime=ifelse(v2x_regime<=1, 0, 1)) %>%
  na.omit %>%
  group_by(year) %>%
  summarize_all(mean, na.rm=T) %>%
  rename(RoW=v2x_regime,
         BMR=e_boix_regime) %>%
  pivot_longer(-year) %>%
  mutate(title="Share of Democratic Regimes") %>%
  ggplot(., aes(x=year, y=value, group=name, color=name, linetype=name)) +
  theme(plot.title=element_text(size=12, hjust=0.5),
        plot.margin=unit(c(0.5, 0.75, 0.5, 0.75), "cm"),
        axis.title.y=element_blank(),
        axis.text.y=element_text(size=11, color="black"),
        axis.title.x=element_blank(),
        axis.text.x=element_text(size=11, color="black"),
        legend.title=element_blank(),
        legend.key=element_rect(colour=NA, fill=NA),
        legend.key.width=unit(1, "cm"),
        legend.text=element_text(size=11),
        legend.position="bottom",
        legend.margin=margin(0,0,0,0),
        legend.box.margin=margin(-5,-5,-5,-5),
        panel.border=element_rect(colour="black", fill=NA, linewidth=1),
        panel.background=element_rect(fill="white"),
        panel.grid.major=element_line(color="grey", linewidth=0.5),
        strip.background=element_rect(color="black", fill=NA),
        strip.text=element_text(size=11, hjust=0.5)) +
  geom_line(linewidth=0.75) +
  scale_linetype_manual(values=c("solid", "dashed")) +
  scale_color_manual(values=c("blue4", "springgreen3")) +
  scale_x_continuous(breaks=seq(1900, 2020, 40),
                     limits=c(1900, 2020),
                     expand=c(0,0)) +
  scale_y_continuous(breaks=seq(0,1,0.25),
                     limits=c(0,1),
                     expand=c(0,0)) +
  facet_grid(.~title)


cp_repression <- vdem %>%
  select(country_id, year, v2x_clphy) %>%
  left_join(., pts %>%
              filter(!(Country=="Czech Republic" & Year<=1991),
                     !(Country=="Czechoslovakia" & Year>=1992),
                     !(Country=="Yugoslavia, Socialist Federal Republic of" & Year>=1993),
                     !(Country=="Yugoslavia, Federal Republic of" & Year<=1992),
                     !(Country=="Yugoslavia, Federal Republic of" & Year>=2007),
                     !(Country=="Serbia" & Year<=2006),
                     !(Country=="Yemen Arab Republic" & Year>=1990),
                     !(Country=="Yemen" & Year<=1989),
                     !(Country=="Soviet Union" & Year>=1992),
                     !(Country=="Russian Federation" & Year<=1991),
                     !(Country=="German Federal Republic" & Year>=1990),
                     !(Country=="Germany" & Year<=1989),
                     !(Country=="Gaza (Hamas)"),
                     !(Country=="Israel in Occupied Territories"),
                     !(Country=="Israel in pre-1967 borders")) %>%
              mutate(Country=ifelse(Country=="Czechoslovakia", "Czech Republic", Country),
                     Country=ifelse(Country=="Yemen Arab Republic", "Yemen", Country),
                     Country=ifelse(Country=="Yugoslavia, Socialist Federal Republic of", "Serbia", Country),
                     Country=ifelse(Country=="Yugoslavia, Federal Republic of", "Serbia", Country)) %>%
              mutate(country_id=countrycode(Country, "country.name", "vdem")) %>%
              filter(!is.na(country_id))  %>%
              rename(year=Year) %>%
              select(country_id, year, PTS_A)) %>%
  left_join(., ciri %>%
              filter(!(CTRY=="Czech Republic" & YEAR<=1992),
                     !(CTRY=="Czechoslovakia" & YEAR>=1993),
                     !(CTRY=="Yugoslavia" & YEAR>=1992),
                     !(CTRY=="Yugoslavia, Federal Republic of" & YEAR<=1999),
                     !(CTRY=="Yugoslavia, Federal Republic of" & YEAR>=2003),
                     !(CTRY=="Serbia and Montenegro" & YEAR<=1991),
                     !(CTRY=="Serbia and Montenegro" & YEAR>=2006),
                     !(CTRY=="Serbia and Montenegro" & YEAR>=2000 & YEAR<=2002),
                     !(CTRY=="Serbia and Montenegro" & YEAR>=2006),
                     !(CTRY=="Serbia" & YEAR<=2005),
                     !(CTRY=="Soviet Union" & YEAR>=1992),
                     !(CTRY=="Russia" & YEAR<=1991),
                     !(CTRY=="Yemen Arab Republic" & YEAR>=1991),
                     !(CTRY=="Yemen" & YEAR<=1990)) %>%
              mutate(CTRY=ifelse(CTRY=="Yugoslavia", "Serbia", CTRY),
                     CTRY=ifelse(CTRY=="Yugoslavia, Federal Republic of", "Serbia", CTRY),
                     CTRY=ifelse(CTRY=="Serbia and Montenegro", "Serbia", CTRY),
                     CTRY=ifelse(CTRY=="Czechoslovakia", "Czech Republic", CTRY),
                     CTRY=ifelse(CTRY=="Soviet Union", "Russia", CTRY),
                     CTRY=ifelse(CTRY=="Yemen Arab Republic", "Yemen", CTRY)) %>%
              mutate(country_id=countrycode(CTRY, "country.name", "vdem")) %>%
              filter(!is.na(country_id)) %>%
              rename(any_of(setNames(names(.), tolower(names(.))))) %>%
              select(country_id, year, physint)) %>%
  mutate(PTS_A=(PTS_A-1)/4,
         physint=physint/8) %>%
  select(-country_id) %>%
  na.omit %>%
  group_by(year) %>%
  summarize_all(mean, na.rm=T) %>%
  rename(PTS=PTS_A,
         CIRI=physint,
         PVI=v2x_clphy) %>%
  pivot_longer(-year) %>%
  mutate(title="Global Mean: Repression") %>%
  ggplot(., aes(x=year, y=value, group=name, color=name, linetype=name)) +
  theme(plot.title=element_text(size=12, hjust=0.5),
        plot.margin=unit(c(0.5, 0.75, 0.5, 0.75), "cm"),
        axis.title.y=element_blank(),
        axis.text.y=element_text(size=11, color="black"),
        axis.title.x=element_blank(),
        axis.text.x=element_text(size=11, color="black"),
        legend.title=element_blank(),
        legend.key=element_rect(colour=NA, fill=NA),
        legend.key.width=unit(1, "cm"),
        legend.text=element_text(size=11),
        legend.position="bottom",
        legend.margin=margin(0,0,0,0),
        legend.box.margin=margin(-5,-5,-5,-5),
        panel.border=element_rect(colour="black", fill=NA, linewidth=1),
        panel.background=element_rect(fill="white"),
        panel.grid.major=element_line(color="grey", linewidth=0.5),
        strip.background=element_rect(color="black", fill=NA),
        strip.text=element_text(size=11, hjust=0.5)) +
  geom_line(linewidth=0.75) +
  scale_linetype_manual(values=c("solid", "dashed", "dotted")) +
  scale_color_manual(values=c("blue4", "springgreen3", "red2")) +
  scale_x_continuous(breaks=seq(1985, 2005, 10),
                     limits=c(1981, 2011),
                     expand=c(0,0)) +
  scale_y_continuous(breaks=seq(0,1,0.25),
                     limits=c(0,1),
                     expand=c(0,0)) +
  facet_grid(.~title)


cp_conflict <-  vdem %>%
  filter(year>=1946) %>%
  select(country_id, year, e_civil_war) %>%
  left_join(., acd %>%
              select(-acd_intra_intensitymax)) %>%
  left_join(., cv %>%
              mutate(country_id=countrycode(CcodeA, "cown", "vdem")) %>%
              mutate(country_id=ifelse(CcodeA==345, 198, country_id),
                     country_id=ifelse(CcodeA==378, 14, country_id)) %>%
              select(-CcodeA)) %>%
  group_by(country_id) %>%
  filter(sum(is.na(acd_intra))<n() & sum(is.na(cv_war))<n()) %>%
  ungroup %>%
  mutate(acd_intra=ifelse(is.na(acd_intra), 0, acd_intra),
         cv_war=ifelse(is.na(cv_war), 0, cv_war)) %>%
  select(-country_id) %>%
  na.omit %>%
  group_by(year) %>%
  summarize_all(mean, na.rm=T) %>%
  rename(HM=e_civil_war,
         UCDP=acd_intra,
         CoW=cv_war) %>%
  pivot_longer(-year) %>%
  mutate(title="Share of Countries with Intra-state Conflict") %>%
  ggplot(., aes(x=year, y=value, group=name, color=name, linetype=name)) +
  theme(plot.title=element_text(size=12, hjust=0.5),
        plot.margin=unit(c(0.5, 0.75, 0.5, 0.75), "cm"),
        axis.title.y=element_blank(),
        axis.text.y=element_text(size=11, color="black"),
        axis.title.x=element_blank(),
        axis.text.x=element_text(size=11, color="black"),
        legend.title=element_blank(),
        legend.key=element_rect(colour=NA, fill=NA),
        legend.key.width=unit(1, "cm"),
        legend.text=element_text(size=11),
        legend.position="bottom",
        legend.margin=margin(0,0,0,0),
        legend.box.margin=margin(-5,-5,-5,-5),
        panel.border=element_rect(colour="black", fill=NA, linewidth=1),
        panel.background=element_rect(fill="white"),
        panel.grid.major=element_line(color="grey", linewidth=0.5),
        strip.background=element_rect(color="black", fill=NA),
        strip.text=element_text(size=11, hjust=0.5)) +
  geom_line(linewidth=0.75) +
  scale_linetype_manual(values=c("solid", "dashed", "dotted")) +
  scale_color_manual(values=c("blue4", "springgreen3", "red2")) +
  scale_x_continuous(breaks=seq(1950, 2010, 20),
                     limits=c(1946, 2006),
                     expand=c(0,0)) +
  scale_y_continuous(breaks=seq(0,1,0.25),
                     limits=c(0,1),
                     expand=c(0,0)) +
  facet_grid(.~title)


fig1 <- ggarrange(cp_conflict, cp_dem, cp_ed, cp_repression, nrow=2, ncol=2)
ggsave(fig1, file="output/Figure 1.pdf", width=8, height=5, dpi=1000)




##############
## Figure 2 ##
##############
fig2a <- main_dat %>%
  select(country_name, year,
         epsm_cent,
         v2xed_ed_cent,
         centralization_HEQ) %>%
  pivot_longer(cols=3:5,
               names_to="var",
               values_to="value") %>%
  mutate(country_name=factor(country_name,
                             levels=c("Spain",
                                      "Italy",
                                      "Germany",
                                      "Chile",
                                      "Argentina"))) %>%
  filter(!is.na(country_name)) %>%
  mutate(var=factor(var,
                    levels=c("epsm_cent",
                             "v2xed_ed_cent",
                             "centralization_HEQ"),
                    labels=c("EPSM", "V-Indoc", "HEQ"))) %>%
  ggplot(., aes(x=year, y=country_name, fill=value)) +
  theme(plot.title=element_text(size=16, hjust=0.5),
        plot.margin=unit(c(0.7, 0.75, 1.625, 0.75), "cm"),
        axis.title.y=element_blank(),
        axis.text.y=element_text(size=14, color="black"),
        axis.title.x=element_blank(),
        axis.text.x=element_text(size=14, color="black"),
        legend.title=element_blank(),
        legend.key=element_rect(colour=NA, fill=NA),
        legend.text=element_text(size=14),
        legend.position="bottom",
        legend.key.width=unit(1, "cm"),
        panel.border=element_rect(colour="black", fill=NA, linewidth=1),
        panel.background=element_rect(fill="white"),
        strip.background=element_rect(color="black", fill=NA),
        strip.text=element_text(size=14)) +
  ggtitle("(A) Centralization Index") +
  geom_tile() +
  scale_fill_gradientn(colors=c("white", "lightblue1", "tomato2"),
                      na.value="black") +
  scale_x_continuous(breaks=seq(1945, 2015, 10),
                     limits=c(1945, 2015),
                     expand=c(0,0)) +
  facet_wrap(~ var, ncol=1)


fig2b <- main_dat %>%
  select(country_name, year,
         edu_power,
         v2edcentcurrlm_ord,
         v2edcentcurrlm_HEQ) %>%
  pivot_longer(cols=3:5,
               names_to="var",
               values_to="value") %>%
  mutate(country_name=factor(country_name,
                             levels=c("Spain",
                                      "Italy",
                                      "Germany",
                                      "Chile",
                                      "Argentina"))) %>%
  filter(!is.na(country_name)) %>%
  mutate(var=factor(var,
                    levels=c("edu_power",
                             "v2edcentcurrlm_ord",
                             "v2edcentcurrlm_HEQ"),
                    labels=c("EPSM", "V-Indoc", "HEQ"))) %>%
  mutate(value=ifelse(is.na(value), -1, value),
         value=factor(value,
                      levels=c(-1, 0, 1, 2),
                      labels=c("NA",
                               "No national authorities",
                               "Both sub-national and national authorities",
                               "National authorities"))) %>%
  ggplot(., aes(x=year, y=country_name, fill=value)) +
  theme(plot.title=element_text(size=16, hjust=0.5),
        plot.margin=unit(c(0.75, 0.75, 0, 0.75), "cm"),
        axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        axis.title.x=element_blank(),
        axis.text.x=element_text(size=14, color="black"),
        legend.title=element_text(size=14),
        legend.key=element_rect(colour=NA, fill=NA),
        legend.text=element_text(size=14),
        legend.position="bottom",
        legend.key.width=unit(1, "cm"),
        panel.border=element_rect(color="black", fill=NA, linewidth=1),
        panel.background=element_rect(fill="white"),
        panel.grid=element_line(color="grey"),
        strip.background=element_rect(color="black", fill=NA),
        strip.text=element_text(size=14)) +
  ggtitle("(B) Centralization of the Curriculum Indicator") +
  geom_tile() +
  scale_fill_manual(breaks=c("No national authorities",
                             "Both sub-national and national authorities",
                             "National authorities"),
                    values=c("white", "lightblue1", "tomato2", "black"),
                    na.value="black") +
  scale_x_continuous(breaks=seq(1945, 2015, 10),
                     limits=c(1945, 2015),
                     expand=c(0,0)) +
  guides(fill=guide_legend(title="Official curriculum set by:",
                           title.position="top",
                           ncol=1)) +
  facet_wrap(~ var, ncol=1)


fig2 <- ggarrange(fig2a, fig2b, nrow=1,  widths=c(1, 0.84))
ggsave(fig2, file="output/Figure 2.pdf", width=10, height=6, dpi=1000)




##############
## Figure 3 ##
##############
fig3 <- main_dat %>%
  filter(country_name %in% c("Argentina", "Chile", "Germany", "Italy", "Spain")) %>%
  select(year, country_name, v2edtehire_ord, polit_train_hire_HEQ) %>%
  gather(var, value, -year, -country_name) %>%
  filter(year<=2015) %>%
  mutate(var=factor(var, levels=c("v2edtehire_ord", "polit_train_hire_HEQ"),
                    labels=c("V-Indoc", "HEQ"))) %>%
  ggplot(aes(x=year, y=value, colour=var, lty=var)) +
  geom_line(linewidth=1.5) +
  scale_color_manual(values=c("blue4", "springgreen3")) +
  facet_wrap(~ country_name, ncol=2) +
  theme(plot.title=element_text(size=11, hjust=0.5, color="black"),
        axis.title.y=element_blank(),
        axis.text.y=element_text(size=11, color="black"),
        axis.title.x=element_blank(),
        axis.text.x=element_text(size=11, color="black"),
        legend.title=element_blank(),
        legend.key=element_rect(colour=NA, fill=NA),
        legend.text=element_text(size=11),
        legend.position="bottom",
        panel.border=element_rect(colour="black", fill=NA, linewidth=1),
        panel.background=element_rect(fill="white"),
        panel.grid=element_line(color="grey"),
        strip.background=element_rect(color="black", fill=NA),
        strip.text=element_text(size=11, color="black")) +
  scale_x_continuous(breaks = seq(1945, 2015, 10), limits = c(1945, 2015)) +
  scale_y_continuous(breaks = seq(0, 1, 1), limits = c(0,1))

ggsave(fig3, file="output/Figure 3.pdf", width=8, height=6, dpi=1000)




##############
## Figure 4 ##
##############
fig4 <- main_dat %>%
  filter(country_name %in% c("Argentina", "Chile", "Germany", "Italy", "Spain")) %>%
  select(year, country_name, national_civics_contents_religion, national_civics_religion_HEQ) %>%
  gather(var, value, -year, -country_name) %>%
  filter(year<=2015) %>%
  mutate(var=factor(var, levels=c("national_civics_contents_religion", "national_civics_religion_HEQ"),
                    labels=c("EPSM", "HEQ"))) %>%
  ggplot(aes(x=year, y=value, colour=var, lty=var)) +
  geom_line(linewidth=1.5) +
  scale_color_manual(values=c("blue4", "springgreen3")) +
  facet_wrap(~ country_name, ncol=2) +
  theme(plot.title=element_text(size=11, hjust=0.5, color="black"),
        axis.title.y=element_blank(),
        axis.text.y=element_text(size=11, color="black"),
        axis.title.x=element_blank(),
        axis.text.x=element_text(size=11, color="black"),
        legend.title=element_blank(),
        legend.key=element_rect(colour=NA, fill=NA),
        legend.text=element_text(size=11),
        legend.position="bottom",
        panel.border=element_rect(colour="black", fill=NA, linewidth=1),
        panel.background=element_rect(fill="white"),
        panel.grid=element_line(color="grey"),
        strip.background=element_rect(color="black", fill=NA),
        strip.text=element_text(size=11, color="black")) +
  scale_x_continuous(breaks = seq(1945, 2015, 10), limits = c(1945, 2015)) +
  scale_y_continuous(breaks = seq(0, 1, 1), limits = c(0,1))

ggsave(fig4, file="output/Figure 4.pdf", width=8, height=6, dpi=1000)




###############
## Figure A1 ##
###############
## Note: data for the figure in the top panel of Figure A1 is not publicly available.

figa1 <- coder_dat %>% 
  group_by(country_id, year) %>% 
  select(ends_with("_conf")) %>%
  select(-v2medstateprint_conf, -v2medstatebroad_conf, -v2medpolstate_conf, 
         -v2medpolnonstate_conf, v2medpatriot_conf, -v2medentrain_conf) %>%
  gather(var, val, -country_id, -year) %>%
  ungroup %>%
  ggplot(., aes(x=year, y=val, group=year)) +
  geom_boxplot(alpha=0.5, color="grey", fill="grey") +
  scale_x_continuous(breaks = seq(1945, 2021, 15), limits = c(1945, 2021)) +
  scale_y_continuous(breaks = seq(0, 1, .25), limits = c(0,1), labels = scales::percent) +
  labs(x="Year", y="Confidence Level, %") +
  theme(axis.title.y=element_text(size=12, color="black"),
        axis.text.y=element_text(size=12, color="black"),
        axis.ticks.y=element_blank(),
        axis.title.x=element_text(size=12, color="black"),
        axis.text.x=element_text(size=12, color="black"),
        legend.position="none",
        legend.title=element_blank(),
        legend.key=element_rect(fill=NA),
        legend.text=element_text(size=12),
        panel.border=element_rect(colour="black", fill=NA, size=1),
        panel.background=element_rect(fill="transparent"),
        panel.grid.major=element_line(color="lightgrey"),
        panel.grid.minor=element_blank(),
        plot.title=element_text(hjust=0.5),
        strip.background=element_rect(color="black", fill=NA),
        strip.text=element_text(size=12, color="black"))

ggsave(figa1, file="output/Figure A1.png", width=7, height=5)    




##################
## Figure B1/B2 ##
##################

## Get data
att_dat <- main_dat %>%
  mutate(cent_index_centralization_eldar_std=scale(epsm_cent),
         cent_index_v2xed_ed_cent_std=scale(v2xed_ed_cent)) %>%
  select(country_id, year,
         edu_power, v2edcentcurrlm_ord,
         cent_index_centralization_eldar_std, cent_index_v2xed_ed_cent_std,
         D2_BMR) %>%
  na.omit



## Run models
dem_var <- c("D2_BMR")
var_list <- c("cent_index_centralization_eldar_std", "cent_index_v2xed_ed_cent_std", "edu_power", "v2edcentcurrlm_ord")


set.seed(2025)
for (i in 1:length(var_list)){
  for (j in 1:length(dem_var)){
    print(paste0(var_list[i]," (", i, "/", length(var_list),")"))
    
    sample_dat <- att_dat %>%
      select(country_id,
             year,
             all_of(var_list[i]),
             all_of(dem_var[j]))
    
    colnames(sample_dat)[3:4] <- c("outcome",
                                   "treat")
    
    out.mc <- fect(outcome~treat,
                   data=sample_dat,
                   index=c("country_id","year"),
                   method="mc",
                   force="two-way",
                   CV=T,
                   na.rm=T,
                   se=T,
                   parallel=T,
                   nboots=1000,
                   nlambda=10)
    
    att.results <- as.data.frame(out.mc$est.att)
    att.results$t <- rownames(att.results)
    write.csv(att.results, paste0("output/results/",var_list[i],"_",dem_var[j],"_att.csv"), row.names=F)
    
    out.plot <- plot(out.mc,
                     main=paste0("ATT: ", dem_var[j], " and ", var_list[i]),
                     stats="F.p",
                     xlim=c(-10,10),
                     cex.main=0.8,
                     cex.lab=0.8,
                     cex.axis=0.8)
    ggsave(paste0("output/results/",var_list[i],"_",dem_var[j],"_att.png"), height=6, width=8)
  }
}



## Figure B1
figb1 <- rbind(
  read.csv(paste0("output/results/",var_list[1],"_D2_BMR_att.csv")) %>%
    filter(t>=-10 & t<=10) %>%
    mutate(period=ifelse(t<0, 0,
                         ifelse(p.value>0.05, 1, 2))) %>%
    mutate(period=factor(period, levels=c(0,1,2), labels=c("Pre", "Post Not Sig", "Post Sig"))) %>%
    mutate(dat_source=1), 
  read.csv(paste0("output/results/",var_list[2],"_D2_BMR_att.csv")) %>%
    filter(t>=-10 & t<=10) %>%
    mutate(period=ifelse(t<0, 0,
                         ifelse(p.value>0.05, 1, 2))) %>%
    mutate(period=factor(period, levels=c(0,1,2), labels=c("Pre", "Post Not Sig", "Post Sig"))) %>%
    mutate(dat_source=2)) %>%
  mutate(dat_source=factor(dat_source, levels=c(1,2), labels=c("Standardized Education Centralization Index (EPSM)",
                                                               "Standardized Education Centralization Index (V-Indoc)"))) %>%
  ggplot(., aes(x=t, y=ATT, color=period)) +
  theme(axis.title.y=element_text(size=12, color="black"),
        axis.text.y=element_text(size=12, color="black"),
        axis.ticks.y=element_blank(),
        axis.title.x=element_text(size=12, color="black"),
        axis.text.x=element_text(size=12, color="black"),
        legend.position="none",
        legend.title=element_blank(),
        legend.key=element_rect(fill=NA),
        legend.text=element_text(size=12),
        panel.border=element_rect(colour="black", fill=NA, size=1),
        panel.background=element_rect(fill="transparent"),
        panel.grid.major=element_line(color="lightgrey"),
        panel.grid.minor=element_blank(),
        plot.title=element_text(hjust=0.5),
        strip.background=element_rect(color="black", fill=NA),
        strip.text=element_text(size=12, color="black")) +
  geom_hline(yintercept=0, color="red") +
  geom_vline(xintercept=0, linetype="dashed") +
  geom_errorbar(aes(ymin=CI.lower, ymax=CI.upper, width=0), size=1) +
  geom_point(size=1.5) +
  scale_color_manual(values=c("black","lightblue", "blue"), drop=F) +
  scale_x_continuous(breaks=seq(-10,10,2)) +
  scale_y_continuous(breaks=round(seq(-0.7, 0.2, 0.1), 1),
                     lim=c(-0.71, 0.21),
                     expand=c(0,0)) +
  xlab("Years since Democratization") +
  ylab("ATT") +
  facet_wrap(vars(dat_source))

ggsave(figb1, file="output/Figure B1.png", width=10, height=4)    
    
    
    
## Figure B2
figb2 <- rbind(
  read.csv(paste0("output/results/",var_list[3],"_D2_BMR_att.csv")) %>%
    filter(t>=-10 & t<=10) %>%
    mutate(period=ifelse(t<0, 0,
                         ifelse(p.value>0.05, 1, 2))) %>%
    mutate(period=factor(period, levels=c(0,1,2), labels=c("Pre", "Post Not Sig", "Post Sig"))) %>%
    mutate(dat_source=1), 
  read.csv(paste0("output/results/",var_list[4],"_D2_BMR_att.csv")) %>%
    filter(t>=-10 & t<=10) %>%
    mutate(period=ifelse(t<0, 0,
                         ifelse(p.value>0.05, 1, 2))) %>%
    mutate(period=factor(period, levels=c(0,1,2), labels=c("Pre", "Post Not Sig", "Post Sig"))) %>%
    mutate(dat_source=2)) %>%
  mutate(dat_source=factor(dat_source, levels=c(1,2), labels=c("Centralization of the Curriculum (EPSM)",
                                                               "Centralization of the Curriculum (V-Indoc)"))) %>%
  ggplot(., aes(x=t, y=ATT, color=period)) +
  theme(axis.title.y=element_text(size=12, color="black"),
        axis.text.y=element_text(size=12, color="black"),
        axis.ticks.y=element_blank(),
        axis.title.x=element_text(size=12, color="black"),
        axis.text.x=element_text(size=12, color="black"),
        legend.position="none",
        legend.title=element_blank(),
        legend.key=element_rect(fill=NA),
        legend.text=element_text(size=12),
        panel.border=element_rect(colour="black", fill=NA, size=1),
        panel.background=element_rect(fill="transparent"),
        panel.grid.major=element_line(color="lightgrey"),
        panel.grid.minor=element_blank(),
        plot.title=element_text(hjust=0.5),
        strip.background=element_rect(color="black", fill=NA),
        strip.text=element_text(size=12, color="black")) +
  geom_hline(yintercept=0, color="red") +
  geom_vline(xintercept=0, linetype="dashed") +
  geom_errorbar(aes(ymin=CI.lower, ymax=CI.upper, width=0), size=1) +
  geom_point(size=1.5) +
  scale_color_manual(values=c("black","lightblue", "blue"), drop=F) +
  scale_x_continuous(breaks=seq(-10,10,2)) +
  scale_y_continuous(breaks=round(seq(-0.3, 0.1, 0.1), 1),
                     lim=c(-0.31, 0.11),
                     expand=c(0,0)) +
  xlab("Years since Democratization") +
  ylab("ATT") +
  facet_wrap(vars(dat_source))

ggsave(figb2, file="output/Figure B2.png", width=10, height=4)




###############
## Figure C1 ##
###############
figc1 <- coder_dat %>%
  filter(country_name %in% c("Argentina", "Chile", "Germany", "Spain", "Italy")) %>%
  distinct(country_name, year, coder_id) %>% 
  group_by(country_name, year) %>%
  summarise(n.coders = n()) %>%
  ggplot(aes(x = year, y = n.coders)) +
  geom_point() +
  facet_wrap(~ country_name)  +
  scale_x_continuous(breaks = seq(1945, 2021, 15), limits = c(1945, 2021)) +
  scale_y_continuous(breaks=seq(0, 10, 3), limits=c(0, 10)) +
  theme(plot.title=element_text(size=11, hjust=0.5, color="black"),
        axis.title.y=element_blank(),
        axis.text.y=element_text(size=11, color="black"),
        axis.title.x=element_blank(),
        axis.text.x=element_text(size=11, color="black"),
        legend.title=element_blank(),
        legend.key=element_rect(colour=NA, fill=NA),
        legend.text=element_text(size=11),
        legend.position="bottom",
        panel.border=element_rect(colour="black", fill=NA, linewidth=1),
        panel.background=element_rect(fill="white"),
        panel.grid=element_line(color="grey"),
        strip.background=element_rect(color="black", fill=NA),
        strip.text=element_text(size=11, color="black"))

ggsave(figc1, file="output/Figure C1.png", width=8, height=4)    





###############
## Figure C2 ##
###############
figc2_plot <- function(plot_var){
  
  p <- coder_dat %>%
    filter(country_name %in% c("Argentina", "Chile", "Germany", "Spain", "Italy")) %>%
    ggplot(aes(x=year, y={{plot_var}}, group=year)) +
    geom_boxplot(alpha=0.5, fill="black", outlier.size=0.5) +
    facet_wrap(~ country_name) +  
    theme(plot.title=element_text(size=11, hjust=0.5, color="black"),
          axis.title.y=element_blank(),
          axis.text.y=element_text(size=11, color="black"),
          axis.title.x=element_blank(),
          axis.text.x=element_text(size=11, color="black"),
          legend.title=element_blank(),
          legend.key=element_rect(colour=NA, fill=NA),
          legend.text=element_text(size=11),
          legend.position="bottom",
          panel.border=element_rect(colour="black", fill=NA, linewidth=1),
          panel.background=element_rect(fill="white"),
          panel.grid=element_line(color="grey"),
          strip.background=element_rect(color="black", fill=NA),
          strip.text=element_text(size=11, color="black")) +
    scale_x_continuous(breaks = seq(1945, 2021, 15), limits = c(1945, 2021)) +
    scale_y_continuous(breaks = seq(0, 1, .25), limits = c(0,1), labels = scales::percent)
  
  return(p)
  
}

figc2a <- figc2_plot(v2edcentcurrlm_conf)
ggsave(figc2a, file="output/Figure C2a.png", width=8, height=6)

figc2b <- figc2_plot(v2edtehire_conf)
ggsave(figc2b, file="output/Figure C2b.png", width=8, height=6)

figc2c <- figc2_plot(v2edideolch_religion_conf)
ggsave(figc2c, file="output/Figure C2c.png", width=8, height=6)













